Digital Twins and Reproducibility

Bridging Point-and-Click COPASI Software with Open-Source and Flexible R Programming language

First Published: Mar 14, 2025

Last updated: mai 19, 2025

Bastien CHASSAGNOL |
Post-Doc
LPSM and Laval

Lilija Wehling |
Ph.D.
VPE BioMed X Institute

Atreyee Banerjee |
Ph.D.
Mainz Universität

Soria Gasparini |
Ph.D.
Heidelberg University

Sanjana Balaji Kuttae |
Computational Biologist (Research Assistant)
VPE BioMed X Institute

Slides available on GitHub

Talks’ Layout

Reproducibility Crisis. 👋
IBD ODE Model 🕘
IBD Numerical Simulations 🕘
BioModels database ⬇️

Reproducibility Crisis

Is there a Reproducibility crisis?

Classes of cell-cell communications cover:

General feeling
Domain-dependant reproducibility
BioModels focus

Nature’s survey of 1,576 researchers on reproducibility in research. Reproduced from Fig1, 1,500 scientists lift the lid on reproducibility, from Baker (2016). — **Nature’s survey of 1,576 researchers on reproducibility in research**. Reproduced from Fig1, *1,500 scientists lift the lid on reproducibility*, from Baker (2016).

Physicists and chemists show the most confidence, in contrast with social or biological studies. Reproduced from Fig2, 1,500 scientists lift the lid on reproducibility, from Baker (2016). — Physicists and chemists show the most confidence, in contrast with social or biological studies. Reproduced from Fig2, *1,500 scientists lift the lid on reproducibility*, from Baker (2016).

Half of public mathematical models stored in BioModels are not curated. Results drawn from 17.05.2025, from Malik-Sheriff et al. (2020). — Half of public mathematical models stored in `BioModels` are not curated. Results drawn from 17.05.2025, from Malik-Sheriff et al. (2020).

\(52\%\) considering there’s a significant crisis
\(90\%\) considering there’s a Reproducibility crisis generally speaking
From Open Science Collaboration (2015) analyses, \(40\%\) of studies in psychology are trustworthy…
… and only \(10\%\) in cancer biology, from Begley and Ellis (2012)

Different types of reproducibility

Empirical reproducibility
Statistical reproducibility
Computational reproducibility

Aka Methodological Reproducibility:

Ability to repeat the same experiment using the same methodology and obtain the same results
Example in scRNA-seq: 2 teams working on the same tissue, get similar distribution of cell type abundances

Reproduce statistical results derived from data analysis (\(p\)-values, effect sizes, confidence intervals, …) using the same statistical methods (or slightly similar ones)
Also named Inferential Reproducibility
Examples in sensitivity analyses: does a change of the random seed or the sampling process induce significant changes of biological findings?¹

Also named Result Reproducibility, can encompass all previously mentioned reproducibility approaches.
Ability to reproduce the computational aspects of a study, which includes reproducing the figures, tables, or other outputs from the data and code provided.
Example in mathematical models: given the same fixed parameters, can you recover the same steady-state conditions?

Reproduce figures and main results from biological papers: Impossible mission?

Main factors contributing to irreproducible research. Unavalaible code and/or datasets is one of the most quoted reasons, while being one of the easiest to address with. Reproduced from Fig.4, Baker (2016).

Mentoring and training emerge out as the main incentives, along with more stringent journal guidelines. Reproduced from Fig.5, Baker (2016).

Repro-hackathons as community-driven events to overcome reproducibility crisis?

Comparison of differentially expressed genes among 12 groups, illustrated by upset plots. Reproduced from Fig.6, Cokelaer, Cohen-Boulakia, and Lemoine (2023) — Comparison of differentially expressed genes among 12 groups, illustrated by **upset plots**. Reproduced from Fig.6, Cokelaer, Cohen-Boulakia, and Lemoine (2023)

Core objective: 12 group students aim at reproducing the main findings from a previously published study, namely the subset of DEGs.
Two-stage Hackathon:
1. Free choice of differential analysis model, genome reference and filtering thresholds (scenario A)
2. Constrained choice of methods and hyperparameters (scenario B)
Conclusion: majority of the differentially expressed genes were common to only a few groups (and no DEG found across all the 12 groups), while all pipelines were runnable.

Repro-hackathon to curate biological models to unified SBML standards

BioModels x BioQuant x BioMedX Team VPE sponsors.

Objective: curate 11 mathematical models of human disease. Cross-disciplinary teams: one biologist, one data scientist, and one computational modeller.
Methods: reproduce at least one output (or recode when only hard-coded ODE equations were provided) of a given biological model + curate species to standard ontologies.
Results: 2 models were completely curated, 4 models were almost recoded and reproduced, 5 models were partially reproduced.

An ODE of the Immune System

Inflammatory Bowel Disease: an ODE model

IBD ODE network is depicted in Figure 1.

Code

flowchart LR
      A["$$x^2$$"] -->|"$$\sqrt{x+3}$$"| B("$$\frac{1}{2}$$")
      A -->|"$$\overbrace{a+b+c}^{\text{note}}$$"| C("$$\pi r^2$$")
      B --> D("$$x = \begin{cases} a &\text{if } b \\ c &\text{if } d \end{cases}$$")
      C --> E("$$x(t)=c_1\begin{bmatrix}-\cos{t}+\sin{t}\\ 2\cos{t} \end{bmatrix}e^{2t}$$")

 flowchart LR
      A["$$x^2$$"] -->|"$$\sqrt{x+3}$$"| B("$$\frac{1}{2}$$")
      A -->|"$$\overbrace{a+b+c}^{\text{note}}$$"| C("$$\pi r^2$$")
      B --> D("$$x = \begin{cases} a &\text{if } b \\ c &\text{if } d \end{cases}$$")
      C --> E("$$x(t)=c_1\begin{bmatrix}-\cos{t}+\sin{t}\\ 2\cos{t} \end{bmatrix}e^{2t}$$")

Figure 1: This is a simple mermaid graph.

Main classes of ODE models, ordered by increasing complexity

First order equations:
Michaelis-Menten equations
Auto-catalytic and Feedback Models
- Positive activation
- negative activation
Complex combination of several processes, pairing together multiple variables

Personalised medecine applied to Inflammatory Bowel Disease

1) Steady-state conditions

Note 1: Derive Steady-States Conditions of an ODE model

Computing the steady-state conditions for an ordinary differential equation (ODE) means finding the system’s equilibrium points, where the variables remain constant over time. Given a system of \(n=15\) ODEs (number of varying species in our model):

\[ \begin{cases} \frac{dx_1}{dt} &= f_1(x_1, x_2, \dots, x_{15}) \\ &\vdots \\ \frac{dx_{15}}{dt} &= f_n(x_1, x_2, \dots, x_{15}) \end{cases} \]

the steady-state conditions are obtained by solving system Equation 1:

\[ \begin{cases} f_1(x_1^*, x_2^*, \dots, x_{15}^*)& = 0 \\ & \vdots \\ f_n(x_1^*, x_2^*, \dots, x_{15}^*) &= 0. \end{cases} \qquad(1)\]

1) Steady-state implementation

healthy_model <- loadModel("../../models/team_2016_final_model_lo2016_2025_05_13.cps")
healthy_model_steady_state <- runSteadyState(
  model = healthy_model,
  calculate_jacobian = TRUE,
  perform_stability_analysis = TRUE,
  method = list("use_newton"=TRUE, accept_negative_concentrations = FALSE)
)

1: Run CoRC::loadModel() function to load the COPASI model in R.
2: We use the CoRC::runSteadyState function to run the steady-states…
3: … and for this task, several algorithms can be chosen, when no explicit solution can be derived directly. We used in our framework the well-known root-finding Newton-Raphson algorithm.

1) Steady-state results

Species	Value(g/cm^3)
M1 macrophage	1.17e-02
M2 macrophage	8.17e-03
Th1 cell	1.40e-01
Th2 cell	3.99e-03
Th17 cell	6.18e-03
Treg cell	2.98e-04
NA	2.01e-07
IL-2	1.07e-08
IL-4	3.36e-09
IL-21	7.77e-08
IL-6	1.46e-06
NA	2.72e-05
IL-10	3.01e-08
NA	1.64e-13
IL-12	3.64e-07
doi:10.1371/journal.pone.0165782.t004

Table 1: Steady-state concentrations of cytokines and immune cells, in healthy individuals in our model.

Original SS reported in Table 4, (Lo et al. 2016, 10). — Original SS reported in **Table 4**, (Lo et al. 2016, 10).

2) Sensitivity analyses

i) Sensitivity analysis Objectives

Sensitivity analysis aims to identify how changes in parameters affect model behaviour, identifying the most influential parameters on the outcome. Besides, these simulation studies contribute to prioritise where uncertainty matters most, and early capture bad representations of the model.

9 parameters were assumed in Lo et al. (2016) to have the strongest impact on the concentrations of T cell subsets. We retrieve their values, inferred from healthy individuals, using the CoRC::getParameters function (see Table 2):

Code

parameters_healthy <- getParameters(model = healthy_model,
                                  key = c("(Diff of M0 to M1).delta_m_cit",
                                          "(Diff of M1 to M2).sigma",
                                          "(Induction of T1 from M).s12",
                                          "(Proliferation of T1).s2",
                                          "(Induction of T2).s4",
                                          "(Induction of T17).s21",
                                          "(Induction of T17).s6",
                                          "(Induction of Tr).sb",
                                          "(Induction of Tr).s10")) |> 
  dplyr::select(-mapping)

flextable(parameters_healthy) |> 
  add_footer_row(values = "key is direct ID of the parameter in the model, reaction describes the biological mechanism associated with the value of the parameter, and value returns the constant value assumed for healthy individuals.", 
                 colwidths = 4) |> 
  bold(part = "header")

key	name	reaction	value
(Diff of M0 to M1).delta_m_cit	delta_m_cit	Diff of M0 to M1	2.40
(Diff of M1 to M2).sigma	sigma	Diff of M1 to M2	24.00
(Induction of T1 from M).s12	s12	Induction of T1 from M	10.93
(Proliferation of T1).s2	s2	Proliferation of T1	1.23
(Induction of T2).s4	s4	Induction of T2	1.94
(Induction of T17).s21	s21	Induction of T17	156.17
(Induction of T17).s6	s6	Induction of T17	156.17
(Induction of Tr).sb	sb	Induction of Tr	14.02
(Induction of Tr).s10	s10	Induction of Tr	14.02
key is direct ID of the parameter in the model, reaction describes the biological mechanism associated with the value of the parameter, and value returns the constant value assumed for healthy individuals.

Table 2: The 9 reaction parameters evaluated for sensitivity analyses.

2) Sensitivity analyses

ii) Latin HyperCube Sampling, step I

Random Sampling for a \(d-\) random vector is performed as Equation 2:

\[ x_i^{(j)} \sim \mathcal{U}\left(x_{\min_j}^{(j)}, x_{\max_j}^{(j)}\right), \quad j = 1, 2, \dots, d, \quad i = 1, 2, \dots, N \qquad(2)\]

Figure 2: **Random sampling**, from Rustell (n.d.). With a low number of observations, the **clustering pattern** inherent to standard random sampling is clearly showcased.

Out of pure randomness, as shown in Figure 2, some regions of the bi-dimensional parameter space remain unsampled, while others are oversampled, as evidenced by visible clustering of points.

In Latin Hypercube Sampling (LHS):

Sampling space is stratified into \(N\) equally-binned intervals for each dimension \(j\)
Permutation is applied across all dimensions to avoid sampling correlation.
Each sample is drawn from a distinct \(N-\) defined \(d-\)sized interval.

Figure 3: **Latin hypercube sample**, from Fig.3 Rustell (n.d.).

LHS provides a better Coverage of the Input Space and Improved Convergence, especially when the sampling size is large with respect to the number of draws, see Figure 3.

Main parameters:

\(x_i^{(j)}\) is the \(i\)-th sample drawn for the \(j\)-th dimension
\(d=9\) (number of parameters included in the sensitivity analysis)
\(N=5000\) the number of parameter distributions simulated
\(\mathcal{U}(x_{\min_j}^{(j)}, x_{\max_j}^{(j)})\) is the Uniform distribution, with bounds defined independently per dimension.

2) Sensitivity analyses

ii) Latin HyperCube Sampling, step II

n_samples <- 5000

parameters_sensitivity <- parameters_healthy$value
param_ranges <- tibble::tibble(parameters_name = parameters_healthy$key,
                               LB= parameters_sensitivity *0.8,
                               UB = parameters_sensitivity *1.2)


set.seed(20)
bootstrap_parameters <- lhs::randomLHS(n_samples, length(parameters_sensitivity))
for (i in seq_along(parameters_sensitivity)) {
  bootstrap_parameters[,i] <- qunif(bootstrap_parameters[,i],
                                    param_ranges$LB[i],
                                    param_ranges$UB[i])
}
apply(bootstrap_parameters,2, quantile)

1: Determine hyper-parameters of sensitivity analysis, such as the number of simulations \(N\), and lower and upper bounds within a \(\pm 20\%\) range.
2: We fix the random seed to ensure reproducibility of the sensitivity analyses.
3: The core function of the lhs package, lhs::randomLHS(), is only able to perform standard uniform sampling.
4: Apply the Inverse Transform Sampling Theorem (see Tip 1) to switch from a standard uniform distribution (\(U \sim \mathcal{U}(0, 1)\)) modelling to any bounded uniform distribution¹
5: Apply the quantile function to verify the quantiles are rather close from their expected theoretical values from uniform sampling.

Tip 1: The Inverse Transform Sampling Theorem

The Inverse Transform Sampling Theorem allows generating any probability distribution from the standard uniform distribution, \(U \sim \mathcal{U}(0,1)\). Let \(X\) be a continuous random variable with cumulative distribution function (CDF) \(F(x)\) and \(F^{-1}\) its reciprocal (the quantile function), then:

\[ X = F^{-1}(U) \]

follows the same distribution as \(X\).

Specifically, by applying the following affine transformation \(X=a + (b - a) U\) on a standard uniformly distributed variable \(u\), \(X\) will follow this distribution: \(\mathcal{U}(a,b)\)

3) Partial Correlation Analyses

Standard Pearson correlation between two variables \(X\) and \(Y\) in Equation 3:

\[ \rho_{X,Y} = \frac{\text{Cov}(X, Y)}{\sigma_X \sigma_Y} \qquad(3)\]

\(\text{Cov}(X, Y) = \mathbb{E}[(X - \mu_X)(Y - \mu_Y)]\) is the empirical covariance matrix.
\(\sigma_X = \sqrt{\text{Var}(X)}\) and \(\sigma_Y = \sqrt{\text{Var}(Y)}\) are standard deviations.

Measures the strength of the linear relationship between \(X\) and \(Y\), without controlling for confounding variables.

The Partial Correlation measures the linear relationship between two variables, conditioned on the remaining set of variables \(\mathrm{Z}\), given by Equation 4:

\[ \rho_{X,Y \mid \mathbf{Z}} = -\frac{\Omega_{X,Y}}{\sqrt{\Omega_{X,X} \Omega_{Y,Y}}} \qquad(4)\]

\(\mathbf{\Omega} = \mathbf{\Sigma}^{-1}\) is the precision matrix.
Non Null off-diagonal terms of the precision matrix reflect direct, causal connections between two variables.

Partial correlation enables more straightforward interpretation by revealing truly causal relationships, while discarding spurious connections.

Practical use case: in Figure 4,

The standard Pearson correlation would have certainly observed a significant correlation between variables \(A\) and \(C\).
This spurious connection between \(A\) and \(C\) due to the confusing effect of \(B\) is clipped using partial correlation instead.

Code

digraph G {
    layout=neato
    A -> B;
    B -> C;
}

Figure 4: Simple example of chain rule associations.

4) Rank Transformation and Scatter Plots

Original Fig.2. reproduced from (Lo et al. 2016, 11)

`BioModels`: a repository of biological models

BioModels: Metadata description

Around half of the publicly avalaible models are tagged **manually-curated**, a ratio slightly increasing over time.

BioModels submission flowchart

From submission to publication of a model in BioModels. Reproduced from Fig4, from Malik-Sheriff et al. (2020)

Conclusion

Take-Aways:

BioModels, a data platform for storing curated biological models.
However, no automated running of code snippets, nor mandatory access to original data
Hard-coded tools, such as COPASI, enable quick analyses for non-expert audience
But they’re limited in terms of functionalities and flexibility, to be used with caution.

This is a Quarto RevealJS Presentation. To learn more, visit https://quarto.org/docs/presentations/revealjs.

If you see mistakes, please raise a GitHub issue.

Bibliographic References

Baker, Monya. 2016. ‘1,500 Scientists Lift the Lid on Reproducibility’. Nature 533 (7604): 452–54. https://doi.org/10.1038/533452a.

Begley, C. Glenn, and Lee M. Ellis. 2012. ‘Raise Standards for Preclinical Cancer Research’. Nature 483 (7391): 531–33. https://doi.org/10.1038/483531a.

Cokelaer, Thomas, Sarah Cohen-Boulakia, and Frédéric Lemoine. 2023. ‘Reprohackathons: Promoting Reproducibility in Bioinformatics Through Training’. Bioinformatics 39 (June): i11–20. https://doi.org/10.1093/bioinformatics/btad227.

Initiative, The Brazilian Reproducibility, Olavo Bohrer Amaral, Clarissa Franca Dias Carneiro, Kleber Neves, Ana Paula Wasilewska Sampaio, Bruna Valerio Gomes, Mariana Boechat de Abreu, et al. 2025. ‘Estimating the Replicability of Brazilian Biomedical Science’. 25 April 2025. https://doi.org/10.1101/2025.04.02.645026.

Lo, Wing-Cheong, Violeta Arsenescu, Razvan I. Arsenescu, and Avner Friedman. 2016. ‘Inflammatory Bowel Disease: How Effective Is TNF-α Suppression?’ PloS One 11 (11): e0165782. https://doi.org/10.1371/journal.pone.0165782.

Malik-Sheriff, Rahuman S, Mihai Glont, Tung V N Nguyen, Krishna Tiwari, Matthew G Roberts, Ashley Xavier, Manh T Vu, et al. 2020. ‘BioModels—15 Years of Sharing Computational Models in Life Science’. Nucleic Acids Research 48 (D1): D407–15. https://doi.org/10.1093/nar/gkz1055.

Open Science Collaboration. 2015. ‘Estimating the Reproducibility of Psychological Science’. Science 349 (6251): aac4716. https://doi.org/10.1126/science.aac4716.

Rustell, Michael. n.d. ‘Knowledge Extraction and the Development of a Decision Support System for the Conceptual Design of Liquefied Natural Gas Terminals Under Risk and Uncertainty.’ University of Surrey. Accessed 14 May 2025. https://openresearch.surrey.ac.uk/esploro/outputs/doctoral/Knowledge-extraction-and-the-development-of/99511011302346.

About

License rights

Opinions expressed are solely my own and do not express the views of the researchers I am associated with.

This work is released under the MIT license.

Made by Quarto Presentations RevealJs

Talks’ Layout

Reproducibility Crisis

An ODE of the Immune System

Personalised medecine applied to Inflammatory Bowel Disease

BioModels: a repository of biological models

Conclusion

Bibliographic References

About

License rights

`BioModels`: a repository of biological models